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Abstract 

A  Monte  Carlo  method  of  solving  the  fundamental  equation  of 
the  kinetic  theory  of  dilute  gases  has  been  developed  and  successfully 
applied  to  two  problems,  one  involving  translational  relaxation  of  a 
spatially  homogeneous  gas  and  the  other  a  plane  steady  shock  of 
arbitrary  strength  (shock  strength  limited  only  by  the  fineness  of  the 
velocity  space  mesh).  This  is  the  first  and  only  method  capable  of 
computing  the  molecular  velocity  distribution  under  conditions  far  from 
equilibrium. 

The  essence  of  the  problem  is  evaluation  of  the  non-linear 
five-dimensional  collision  integral.  Straight  forward  numerical 
quadrature  would  require  about  a  year  on  the  fastest  present  day 
computers.  The  computation  time  is  reduced  to  a  practical  value,  of 
the  order  of  an  hour,  by  a  statistical  sampling  technique  closely 
resembling  the  real  statistical  collision  phenomena  in  the  gas. 

Computations  to  date  have  been  restricted  to  elastic  sphere 
scattering  of  molecules  without  internal  degrees  of  freedom.  Differ¬ 
ential  scattering  cross-sections  other  than  elastic  sphere  can  be 
accommodated  in  the  computer  program  without  complications  or  computing 
time  penalty.  Introduction  of  one  or  two  internal  molecular  degrees 
of  freedom  will  increase  the  complexity  and  computing  time,  but  not 
to  an  impractical  degree. 

Several  technical  problems  had  to  be  solved  in  order  to 
make  the  method  work  properly.  The  "fairness"  of  the  random 
collision-selecting  process  had  to  be  designed  with  care  and  verifie1 
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thoroughly.  Furthermore,  a  tendency  for  the  theoretically  conserved 
quantities  (number  of  molecules,  momentum  and  energy)  to  change 
slightly  because  of  interpolation  and  quadrature  errors,  had  to  be 
corrected  to  prevent  its  building  up  significantly  over  many  time 
steps  or  iterations.  In  the  case  of  the  mere  difficult  shock  wave 
computation,  a  tendency  for  the  shock  front  to  creep  out  of  the 
computational  reference  frame  in  the  course  of  successive  iterations 
had  to  be  eliminated  by  choosing  density  rather  than  distance  normal 
to  the  shock  front  as  independent  variable.  Finally,  a  verifiably 
convergent  iteration  scheme  had  to  be  devised  for  successively  improving 
an  initial  trial  solution.  All  of  these  technical  problems  have  been 
solved.  The  initial  trial  solution  so  far  used  has  been  the  Mott- 
Smith  approximation. 

We  exhibit  for  the  translational  relaxation  problem  a  graph 
of  the  temporal  behavior  of  the  Boltzmann  function  for  Mach  numbers 
ranging  from  0.5  to  6.  For  the  shock  wave  problem  we  show,  for  Mach 
number  3.0,  contour  plots  of  the  Mott-Smith  velocity  distribution 
function,  and  of  the  collision  integral  derived  from  it,  together 
with  other  plots  characterizing  the  Monte  Carlo  solution  of  the 


problem. 
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1.  Introduction 


The  basic  equation  in  the  kinetic  theory  of  dilute  gases  is 
the  Boltzmann  equation.  Analytical  solution  of  this  equation  has  only 
been  possible  for  conditions  near  equilibrium  and  for  a  limited 
number  of  intermolecular  potentials. 

The  principal  obstacle  in  the  way  of  solving  the  Boltzmann 
equation  is  the  nonlinearity  and  complexity  °f  the  collision  integral 
in  it.  Evaluation  of  the  integral  by  direct  numerical  quadrature 
is  too  slow  to  be  useful,  even  on  the  fastest  computers.  The  senior 
author  therefore  devised  in  1955  a  Monte  Carlo  method*-  for  the  numerical 
evaluation  of  the  integral  on  digital  computers.  The  method  can 
be  applied  to  any  velocity  distribution  function.  We  have  programmed 
the  method  for  elastic  sphere  molecules  but  extensions  to  other  molecular 
forces  are  feasible  whenever  the  differential  cross-sections  are  known. 

Since  1955,  the  authors  have  refined  and  carefully  tested 

* 

this  method  and  have  applied  it  to  several  problems  in  kinetic  theory. 

The  present  paper  gives  the  first  complete  account  of  the  method.  We 
made  careful  studies  of  accuracy  which  are  especially  important  for  our 
work  because  of  the  complexity  of  the  collision  integral  and  the  wide 
usefulness  of  a  reliable  means  of  evaluating  it.  All  indications  are 
that  the  Monte  Carlo  calculations  are  valid  to  within  the  expected 


These  studies  were  first  described  in  nine  CSL  reports  published 
in  the  years  1962-1966. 
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statistical  fluctuations.  Lack  of  space  prohibits  detailed  description 
of  the  accuracy  studies  here. 

One  problem  in  kinetic  theory  to  which  we  have  applied  our 

Monte  Carlo  method  is  that  of  the  "pseudo-shock,"  a  kind  of  translational 

2  3 

relaxation  of  molecules.  ’  Solution  of  the  Boltzmann  equation  follows 

easily  once  the  Monte  Carlo  method  permits  evaluation  of  the  necessary 

collision  integrals,  and  the  solution  is  obtained  with  equal  ease  for 

conditions  very  near  and  very  far  from  equilibrium. 

In  a  second  problem,  that  of  shock  structure,  a  more  difficult 

iteration  procedure  is  necessary  to  solve  the  Boltzmann  equation. 

We  have  substantial  numerical  evidence  of  the  convergence  of  the 

iteration  process.  Even  apart  from  the  convergence  question,  the 

Monte  Carlo  evaluation  of  the  collision  integral  can  yield  new  and 

fundamental  information  about  shock  structure.  In  particular  it 
4 

is  now  possible;  a)  to  test  any  velocity  distribution  function 

proposed  as  an  approximate  solution  of  the  Boltzmann  equation  for  a 

shock  wave  or  other  /low  condition;  and  b)  to  check  directly  the 

various  elaborate  analytical  calculations  involved  in  moment  methods. 

5  6  7  8 

The  Monte  Carlo  methods  of  Haviland  ’  and  Bird  ’  have 
also  been  used  in  kinetic  theory  problems.  Neither  method  evaluates 
collision  integrals,  so  neither  method  can  yield  detailed,  accurate, 
and  explicit  solutions  of  the  Boltzmann  equation.  Bird's  method  is, 
however,  nicely  complementary  to  ours  in  that,  though  less  accurate, 
it  is  at  present  fast°r  than  ours  ap-J  is  therefore  already  useful  in 
problems  involving  more  than  one  independent  variable. 
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2.  Outline  of  the  Method 

In  the  Monte  Carlo  method  the  collision  integral  is  first 

replaced  by  an  integral  over  a  finite  region  of  velocity  space.  This 

or  similar  approximations  must  always  be  used  in  numerical  evaluation 

of  integrals  over  an  infinite  region.  The  finite  region,  of  volume  R, 

is  taken  large  enough  so  that  it  includes  most  of  the  molecules.  The 

average  of  the  integrand  over  R  and  over  all  values  of  the  line-of- 

centers  vector  ic  is  then  approximated  by  the  average  of  a  large  and 

fair  sample  of  N  values  of  the  integrand.  A  Monte  Carlo  estimate  of 

the  value  of  the  collision  integral  (with  random  errors  proportional 
-1/2 

to  N  )  is  given  by  the  product  of  this  average  value  with  the 
volume  R.  Note  that  the  integrand  is  a  function  of  eight  independent 
variables  derived  from  the  three  vectors,  v,  v* ,  and  ic  that  define  a 
collision.  Nordsieck's  Monte  Carlo  method  enforces  a  fairness  of  the 
sampling  in  the  eight-dimensional  space  of  these  variables. 

3.  The  Boltzmann  Equation 

The  units  we  shall  use  are  the  values,  denoted  by  the 
subscript  1,  of  various  properties  of  a  reference  gas.  Thus  n^ , 

T^  are  the  units  of  number  density  n  and  temperature  t.  The  unit 
length  of  length  i =  l/(2rm^a  )  =  (mean  free  path)^/  */2.  The  unit 
of  velocity  c^  =  ,/(2nkT^/m)  =  (mean  speed)^  X  (tt/2).  The  unit  of 
time  is  therefore  (mean  free  time)^  X  (/2/rr )  and  of  the  velocity 
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distribution  function  is  n^/c^  • 


In  these  units  the  Boltzmann  equation 
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may  be  written 

df/dT  +  v^df/dx  =  a  -  £f  = 

/(FF'-ff’)  |k*vr|dv'(dk/4TT>.  (1) 

where  f  =  f(v,x,T)  is  the  velocity  distribution  function,  x  is  the 
distance  variable,  and  t  is  the  time  variable;  the  unit  vector  k  gives 
the  direction  of  the  line  of  centers  during  a  collision;  v^  =  v*  -  v; 
and  f ,f’ ,F,F'  denote  the  four  values  of  f  corresponding  to  the  four 
velocities,  v,v',V,V'.  Integration  is  over  the  whole  4rr  solid  angle  in 
order  that  the  k  integration  limits  may  be  independent  of  v  and  v*.  The 
notation  £f  reminds  us  that  this  second  part  of  the  collision  integral 
is  proportional  to  f(v,x,T),  a  fact  of  importance  in  devising  a  stable 
method  of  integrating  the  differential  equation  for  the  case  of  the 
strong  shock  wave. 

4.  Geometrical  and  Numerical  Assumptions 

To  develop  a  specific  Monte  Carlo  algorithm  that  will  yield 
estimates  of  the  Boltzmann  collision  integral  we  must  make  a  number 
of  geometrical  and  numerical  assumptions.  Most  of  these  may  be  easily 
modified  as  required  in  the  future.  First,  as  we  restrict  ourselves  to 
flows  possessing  axial  symmetry,  we  represent  any  velocity  vector  v 
by  its  cylindrical  components  (v  ,v  ,0)  in  velocity  space.  Because 

X  J_ 

of  the  axial  symmetry  the  velocity  distribution  function  will  not 
depend  upon  0.  In  the  sampling,  however,  we  must  still  treat  a  collision 
as  a  full  three-dimensional  phenomenon. 
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Convenience  of  computer  progressing  of  the  Monte  Carlo  scheme 
requires  that  the  region  R  in  velocity  space  be  spanned  by  a  fixed  set 

of  velocity  cells.  We  therefore  transfora  v  to  a  new  variable  v  bv 

o' 

introducing  as  adjustable  parameters,  a  shift  of  the  origin  and  a 
velocity  scale  factor  For  any  given  problem  these  parameters  are 

fixed  so  that  all  but  0.1%,  for  example,  of  the  oolecules  are  within  R. 
Because  of  the  rapid  (Gaussian)  decrease  of  density  with  increase  of 
velocity,  quadrature  error  in  the  velocity  space  is  quite  insensitive 
to  the  fraction  of  oolecules  excluded  from  R  so  long  as  that  fraction 
is  small . 

The  quantization  of  the  velocity  space  v  in  which  we  define 

□ 

fixed  cells  was  designed  to  yield  accuracy  of  the  order  of  1%  in  the 

Monte  Carlo  estimation  of  the  two  parts  of  the  collision  integral. 

(The  values  of  the  components  of  v  and  k  that  are  used  in  the  caicu- 

m  m 

lations  are  those  corresponding  to  the  centers  of  the  cells  in  v  space 

m 

and  space.)  We  choose  226  cells  in  the  two-dimensional  (VX»^L)  space, 

as  shown  in  Fig.  1,  to  cover  the  semicircular  region  v  <  24  in  such 

m 

a  way  that  the  area  of  the  226  cells  is  different  from  the  area  of  the 
semicircular  region  by  less  than  0.1%.  In  our  calculations  we  are 
concerned  with  those  functions  a,  bf,  and  f  that  depend  only  upon  the 
two  velocity  components  v  and  v  ,  together  with  either  the  time  or 
position  variables.  Therefore,  tables  of  226  values  of  each  of  these 
functions  are  stored  in  the  computer  for  each  value  of  the  time  or 
position  variable.  As  examples  of  the  graphical  representation  of 
such  functions  we  show,  in  Fig.  2.  isolines  of  the  Mott-Smith  approximate 
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velocity  distribution  and  the  associated  collision  integral  within  a 
strong  shock  wave.  The  v^  axis  in  the  figures  is  perpendicular  tc 
the  shock  front. 

The  quantized  values  of  the  azimuthal  angle  0  are  chosen  to 

be  odd  multiples  of  (90°/32)  so  that  the  range  of  0  from  0  to  2n 

contains  64  cells.  Sixteen  values  each  of  sin  0  and  cos  0  are  stored 

in  fixed  tables  corresponding  to  the  range  0  to  tt/2.  The  ether 

quadrants  are  introduced  in  the  Monte  Carlo  calculations  by  randomizing 

signs  appropriately.  The  unit  vector  k  is  represented  by  64  sets  of 

values  of  the  three  direction  cosines.  The  values  are  chosen  in 

such  a  way  as  to  divide  the  first  octant  of  the  unit  sphere  into 

64  equal  parts.  These  192  values  of  the  direction  cosines  appear 

in  a  fixed  table  in  the  computer. 

Having  defined  the  cells  we  shall  use  in  v  and  Ic  space,  we 

can  now  define  the  sampling  algorithm.  A  sequence  of  pseudo-random 

,  10 

numbers  each  containing  37  bits  is  generated  by  a  modified  Juncosa 

process.  Bits  in  the  same  number  and  in  the  sequence  are  practically 

37 

uncorrelated  over  the  whole  repeat  cycle  of  2  numbers.  From  each 
random  number  we  derive  a  random  collision  by  using  the  bits  as 
follows:  14  bits  for  the  vector  v,  14  for  the  vector  v' ,  and  9  for 
the  unit  vector  k.  Successive  numbers  in  the  sequence  then  lead  to 
successive  and  independent  collisions  randomly  chosen  from  the 
’  (226/256)^  _  ^  x  cells  in  the  eight  dimensional 

(v,v',k)  space.  A  collision  is  rejected  as  "unsuccessful"  if  either 
V  or  V'  ,  computed  from  (v^'jk),  lies  outside  the  region  Vm  <  24. 
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The  fraction  of  collisions  thus  rejected  is  0.1583.  The  val«-  s  of 

(V  V  ,V  ' ,  V  ')  are,  in  effect,  rounded  to  values  corresponding  to 
x  J-  x  m 

the  nearest  cell  center  in  (v  ,v  )  space. 

X  JL.  Si 


5.  Various  Monte  Carlo  Estimates  of  the  Collision  Integral 


5.1  The  Eight  Basic  Estimates.  Let  us  now  consider  four 
Monte  Carlo  approximations  to  each  of  the  functions  av^  and  £fv 
The  eight  approximating  functions  are: 


a  =  J^FF'h.)* 
~  I  1  . 


bf  =  J  <££'k  >*v 
v 


a*  =  J  (FF'h  )aV 
1  1  v' 


b’f'  =  J.  (f f  * H-  )aV 
~  1  1 

v* 


A  =■  J  <££'*  >^V 
V 


BF  *  J,<FF'x  )3' 

1  I  " 

*  1  tT 


A'  =  J  (ff ’k  >3V 

#«S-r  1  I  — 

L  v! 


B'F'  =  J  (FF* x, > 

1  1  v' 


In  these  equations 


=■■  (8n)  226(K1)"4;  ^  =  [v±v±'  |fc*(v'-v)  |]t 


where  we  are  using  the  "machine  units"  defined  above.  The  symbol  below 
each  average-value  sign  in  Eq.  2  indicates  which  of  the  four  velocities 
is  held  constant  in  averaging  over  a  sample  of  collisions.  Each  of 
the  functions  a, a1, A, A'  approaches  the  function  a,  and  each  of  the 
functions  bf,  b’f',  BF,  B'F'  approaches  the  function  {>f  as  one  decreases 


■rr~c  *  **«*~*vt 
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by  one  means  or  another  the  errors  of  various  kinds.  For  finite  errors 
of  each  type,  however,  the  four  functions  of  each  type  are  distinct. 

We  may  expect  that  the  statistical  errors  in  the  Monte  Carlo 
method  of  evaluating  a  and  ftf  will  decrease  in  proportion  to  N 
where  N  is  the  sample  size.  We  may  expect  that  the  quadrature  error, 
corresponding  to  the  finite  size  of  the  cells  in  v  and  k  spaces  will 
decrease  as  the  volumes  of  the  cells  are  decreased.  Our  program  permits 
easy  variation  of  the  cell  size  for  the  velocity  space.  The  error 
owing  to  the  neglect  of  a  few  molecules  with  large  velocities  will 
decrease  rapidly  as  is  decreased. 

5.2  Linear  Combinations  of  the  Estimates.  An  important 
characteristic  of  the  Monte  Carlo  method  is  that  all  of  the  eight 
functions  in  Eq.  2  can  be  evaluated  simultaneously  and,  indeed,  in 
any  useful  linear  combination.  For  example,  let  us  look  at  the 
evaluation  of  a  which  is  found  by  forming  a  sum  over  the  collision 
sample 

<FF'K.)aV  «  M_1  FF'k  (4) 

1  -  o -  1 

V  V 

where  v  is  held  constant  during  the  summation  and  M  is  equal  to  the 

o 

number  of  increments  or  hits  for  each  velocity  cell  during  the  sampling 
process.  Each  such  sum  is  formed  by  fixed  point  arithmetic  to  permit 
unbiased  rounding  of  the  increments. 

An  increment  (FF'h^/M^)  in  the  value  of  the  sum  is  made, 
for  each  random  sampled  collision  (v ,v' ,k)-*(V, V' ) ,  to  the  velocity 
cell  labelled  by  the  velocity  v.  In  forming  bf  or  b'f'  we  sum 


f 

f 
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(ff'K^/M  )  similarly,  again  making  increments  in  the  v  or  v*  cell, 

respectively.  In  forming  A,  A',  BF,  B’F',  we  make  use  of  the  fact 

^ 

that  each  collision  (for  molecules  exerting  central  forces)  has  an 
inverse;  that  is,  if  (v,v'  >k)-*(V,V’)  is  a  collision,  so  also  is 
(V,V*  ,-ic)-*(v,v')  •  Then  to  get  A  we  sum  (ff’K./H  )  over  (v,v'  k) 

<v  i  o 


while  holding  constant,  and  similarly  for  A',  BF,  and  B'F' . 

In  the  simultaneous  calculation  of  the  eight  quantities  in 

Eq.  2  we  symmetrize  the  Monte  Carlo  summand  by  multiplication  by 

v^.  This  accomplishes  two  ends;  it  increases  the  speed  of  calculating 

the  eight  estimates  and  it  makes  possible  achievement  of  good  statistics 

for  the  cells  near  v^  =  0  because  we  can  distribute  hits  uniformly 

over  (v  ,v.  )  space  and  then  later  correct  by  the  factor  2ttv’  to 
x  JL  4- 

effect  uniform  a  priori  probability  of  collision  over  three  dimensional 
velocity  space. 

The  eight  functions  defined  in  Eq.  2  are  basic  to  all 
considerations  of  Nordsieck’s  Monte  Carlo  method.  In  the  calculations 
discussed  iu  this  paper,  and  in  our  study  of  accuracy,  we  considered 
in  particular  three  linear  combinations  of  these  basic  Monte  Carlo 
functions.  Two  of  these  are 


a  =  l/2(a  +  a’) ,  (5) 

a  =  l/2[a  (v  ,v  )  +  a  (-v  ,v  )],  (6) 

with  similar  equations  for  ,  (bf)^,  ...  ,  BF.  Eq .  5  expresses  an 
average  over  the  two  incoming  (or  two  outgoing)  molecules.  Eq .  6 
expresses  an  average  of  a  Monte  Carlo  estimate  of  a  function  and  of  the 


"reflected"  estimate.  In  Nordsieck's  Monte  Carlo  method  the  samples 
of  v  and  cf  v*  are  ea  .h  distributed  uniformly  and  independently  over 
the  (v  ,v,  )  space  (v  <  24).  In  our  studv  of  hits  per  cell  we  then 
needed  to  consider  only  a  ,  and  not  a  and  a*  separately. 

U  ^  rs* 

The  other  linear  combination  used  is 

aAN  "  1/2<3  +  A) 

with  a  similar  equation  for  (bf) ,v .  This  is  an  average  of  the  Monte 

AN 

Carlo  functions  for  the  direct  and  inverse  collisions  which  takes 
advantage  of  the  uniformity  of  sampling  in  v  and  V  spaces. 

Each  of  the  linear  combinations  that  has  been  defined  is 
relevant  to  the  Monte  Carlo  method  because  the  corresponding 
summands  in  each  case  can  be  computed  from  the  two  basic  functions 

(ff'H  )  and  (FF’h  )  with  almost  no  additional  arithmetic.  Therefore . 

JL  1  » 

the  size  of  a  Monte  Carlo  sample  can  be  increased,  in  effect,  by 
these  exchange  operations  between  primed  and  unprimed  molecules. 
between  incoming  and  outgoing  molecules,  and  between  a  collision  and 
its  reflection. 

5.3  Corrections  that  Force  Conservation.  The  starting 
noint  in  deriving  these  corrections  is  the  observation  that  accurate 
solutions  of  the  Boltzmann  equation  must  satisfy  conservation 
equations  characteristic  of  the  problem  at  hand.  We  therefore  compute  the 
smallest  corrections  of  the  values  of  the  226  values  of  the  collision 
integral,  in  a  least-squares  sense,  that  are  consistent  with  enforcing 


conservation. 
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Why  is  this  correction  method  appropriate?  Ar.y  numerical 
'quadrature  method  of  evaluating  the  collision  integral  produces  non¬ 
zero  quadrature  errors.  There  will  consequently  be  a  bias  in  the 
computed  values  of  the  two  parts  of  the  collision  integral.  (The 
neglect  of  a  few  molecules  with  large  speeds  will  cause  a  somewhat 
smaller  Mite.)  Solutions  of  the  Boltzmann  equation  both  for  the  pseudo¬ 
shock  and  the  shock  wave  require  repeated  evaluation  of  the  collision 
integral,  so  that  uncorrected  bias  tends  to  build  up  unacceptable 
trends  in  the  computed  solutions.  If,  as  in  the  case  of  our  Monte 
Carlo  method,  the  bias  is  not  large,  then  the  least  square  adjustment 
that  we  have  described  will  produce  values  of  the  collision  integral 
which  satisfy  the  conservation  equations  and  prevent,  to  a  large 
extent,  the  drift  or  trend  in  the  calculations. 

6.  The  Computer  Program 

As  noted  above,  our  program  at  present  assumes  elastic 
sphere  molecules  and  axial  symmetry  of  velocity  distributions.  The 
first  restriction  is  easily  removable.  The  program  can  be  modified 
to  remove  the  second  restriction.  However,  the  present  program 
without  modification  is  directly  applicable  to  the  pseudo-shock,  the 
shock  wave  and  all  similar  problems  involving  only  one  space  or  time 
variable . 

The  Monte  Carle  part  of  the  program  contains  about  900  words 
(order-pairs,  constants,  and  fixed  tables).  The  complete  program 
contains  about  8,000  words  which  includes  the  Monte  Carlo  part  of  the 
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program,  integration  and  moment  programs,  and  many  useful  auxiliary 

programs-  As  noted  earlier,  the  number  of  cells  in  the  eight  dimensional 

(v,v’,k)  space  that  is  sampled  is  more  than  10^^.  (A  direct  simple, 

9 

numerical  quadrature  of  comparable  accuracy  would  require  ^  x  10  cells 

in  a  seven  dimensional  space.)  The  program  computes  226  values  (in 

parallel)  of  each  of  the  functions  and  (bf).,,,  according  to  the 

AN  AN 

Nordsieck  prescription,  in  three  minutes  on  the  CDC  1604  with  statistical 
fluctuations  of  8%. 

In  any  accurate  method  of  solution  of  the  Boltzmann  equation 
it  is  necessary  to  add  efficient  methods  of  monitoring  and  judging 
the  large  volume  of  numerical  data  produced.  Our  calculations  generally 
produce  226  values  of  each  of  the  functions  f,  a,  and  bf  for  each 
iteration  in  the  shock  wave  problem  or  each  step  in  time  in  a  relaxation 
problem  like  the  pseudo-shock. 

We  have  two  monitoring  devices  for  following  this  large 

■k 

output  of  data.  One  monitoring  device  is  the  calculation  and  output 
of  many  moments  of  f,  a,  or  (a-bf) .  Our  numerical  integrations  in 
velocity  space  are  good  to  much  better  than  1%,  except  for  high  order 
moments  or  sharply  peaked  integrands.  We  output  11  independent  moments 
of  the  velocity  distribution  function  f  (including  the  two  or  three 
moments  which  should  be  conserved,  the  density,  the  Boltzmann  function, 


It  must  be  remembered  that  in  developing  a  method  of  solving  a  difficult 
equation  like  the  Boltzmann  equation  it  is  necessary  to  make  many  more 
trials  and  test  calculations  than  are  needed  for  the  final  production 
run. 


13 


and  its  flux).  We  also  output  nine  moments  of  a  and  of  (a-bf),  some 
of  which ,  by  conservation,  should  be  zero,  or,  by  the  Boltzmann  theorem, 
should  be  negative. 

As  a  second  monitoring  device  we  output  characters  denoting 
ranges  of  values  of  each  of  the  velocity  dependent  functions  f,  a, 
and  bf.  Approximate  isolines  may  easily  be  derived  from  these  characters 
and  permit  immediate  and  vivid  visualization  of  the  nature  of  the 
functions  and  their  changes  from  one  physical  or  calculational  situation 
to  another.  The  spacing  of  the  isolines  can  be  set  at  any  appropriate 
level  by  changing  the  table  of  the  reference  function  values  that 
correspond  to  the  characters  output.  Examples  of  accurate  isolines, 
derived  by  such  a  procedure,  are  shown  in  Fig.  2. 

7.  Behavior  of  Solutions  of  the  Boltzmann  Equation 

We  have  been  particularly  interested  in  using  the  Monte 
Carlo  evaluation  of  the  Boltzmann  collision  integral  in  calculating 
a  translational  relaxation  process  (the  "pseudo-shock")  and  in  trj’ing 
to  solve  the  Boltzmann  equation  for  a  shock  wave.  From  each  of  these 
problems  we  have  derived  considerable,  though  indirect,  support  for 
the  validity  and  reliability  of  the  Monte  Carlo  method.  In  both 
problems  we  takf-,  as  absolute  requirements  on  the  validity  of  the 
method,  that  the  appropriate  conservation  equations  be  satisfied 
and  that  the  appropriate  Boltzmann  function  decrease  monotonically . 

Such  requirements  must  clearly  be  made  on  physical  grounds  and  are 
especially  important  in  the  absence  of  methods  other  than  Monte 
Carlo  of  solving  the  Boltzmann  equation  for  these  problems. 
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Using  the  least  square  adjustments  of  the  collision  integrals 

3 

we  solved  the  pseudo-shock  problem  and  demonstrated  the  accuracy, 

stability,  and  convergence  of  the  solution  process,  these  characteristics 

being  dependent  upon  the  reliability  of  the  Monte  Carlo  evaluation  of 

the  collision  integral  Fig.  3  illustrates  several  of  the  favorable 

characteristics  of  our  Monte  Carlo  solution  for  this  pseudo-shock. 

The  Boltzmann  function  decreases  monotonicaliy ,  for  each  value  of  the 

Mach  number  M,  and  Monte  Carlo  fluctuations  have  negligible  effect 

until  late  in  the  relaxation  process.  Even  for  the  rather  small  samples 

used  (10,000  collisions  for  each  forward  step  in  time)  we  were  able 

to  calculate  the  slope  of  each  of  the  curves  in  Fig.  3  and  show  that 

there  was  no  significant  difference  between  the  slopes. 

It  has  been  our  objective  since  1958  to  solve  the  Bol -zmann 

equation  for  shock  waves  using  our  Monte  Carlo  method.  It  has  been 

apparent  for  several  years  that  it  is  the  convergence  and  stability  of 

our  iterative  method  of  solution  of  the  Boltzmann  differential  equation, 

rather  than  the  difficulty  in  evaluating  the  collision  integral  in  it, 

which  has  been  impeding  our  progress  toward  solution  of  this  problem. 

We  have  recently  been  using  the  correction  method  which  forces  the 

conservation  equations  to  hold  and  have  begun  obtaining  stable  and 

convergent  results  in  treating  the  internal  structure  of  shock  waves. 

Let  us  illustrate  these  remarks  with  data  for  a  shock  wave 

that  is  described  by  the  values  of  f(v,x)  and  of  (a-bf)  at  nine 

positions  within  the  shock.  Fig.  4  illustrates  the  rapid  convergence 

13 

obtained  for  a  typical  Monte  Carlo  run  for  M  =  2.5  with  2  collisions 


2 


in  each  Monte  Carlo  sample  used.  The  rtns  change,  5f,  of  f  from  one 
iteration  to  the  next,  decreases  by  a  factor  of  seven  in  the  first 
three  iterations.  By  the  twelfth  iteration  the  rms  value  of  6f  has 
decreased  to  0.00013,  although  the  largest  value  of  f,  near  the  cold 
boundary  of  the  shock,  is  close  to  unity. 

The  data  shown  in  Fig.  5  were  obtained  by  averaging  the 
values  from  four  independent  runs  of  twelve  iterations  each,  with 


2  collisions  in  each  Monte  Carlo  sample.  It  has  been  shown 
earlier^  that  the  density  profile  n’  is  appropriately  represented 
by  a  plot  of  the  density  gradient  dn/dx  against  the  reduced  number 


density  n^n-n^/fa^-n^) ,  as  in  Fig.  5.  For  a  Mott-Smith  shock 

the  curve  is  a  parabola.  Even  at  this  exploratory  stage  of  our 
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calculations,  with  samples  of  2  ,  there  is  some  evidence  that  the 

shock  is  significantly  asymmetrical.  As  expected  from  the  Boltzmann 
theorem,  the  Boltzmann  flux  derived  from  the  same  four  runs  decreases 
monotonically  through  the  shock. 

Further  indirect  checks  of  the  Monte  Carlo  evaluation  of  the 
collision  integral  will  depend  upon  comparisons  with  accurate 
analytical  solutions  of  the  Boltzmann  equation  for  very  weak  or  very 
strong  shockwaves.  The  Navier-Stokes  equations  should  describe 
accurately  the  internal  structure  of  weak  shock  waves,  and  also  the 

characteristics  of  shocks  of  any  strength  near  their  up-  and  down-stream 

11  12 
boundaries.  Although  it  has  been  suggested  that  the  Mott-Smith 
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model  is  asymptotically  correct  for  strong  shocks,  it  predicts  the 
wrong  value  of  the  Prandtl  number  near  the  down-stream  boundary. ^ 
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In  view  of  the  errors  of  the  various  methods  we  do  not  know  at  this 
time  whether  there  is  any  value  of  the  Mach  number  for  which  the  errors 
of  the  methods  would  be  small  enough  so  that  comparisons  of  the  predicted 
details  of  shock  structure  would  be  meaningful.  Efforts  to  make 
such  meaningful  comparisons  are  obviously  important  both  to  check 
further  the  solutions  based  on  the  Monte  Carlo  evaluation  of  the  collision 
integral  and  to  define  the  range  of  validity  of  the  analytical  theories. 

8.  Future  Studies 

We  shall  outline  here  briefly  a  variety  of  directions  that 
may  be  taken  in  future  research,  now  that  there  is  a  reliable 
method  of  evaluating  the  Boltzmann  collision  integral.  These  various 
researches  may  involve  either  modifications  of  our  method  or 
applications  of  it,  and  these  aspects  of  course  will  often  overlap. 

Let  us  first  consider  modifications  and  studies  of  the 
Monte  Carlo  method  which  will  be  useful  in  all  of  the  applications  of 
this  method  that  will  be  outlined  below.  These  are:  a)  Many  random 
collisions  calculated  once  for  all  and  prestored  on  tape  or  disc. 

Use  of  such  prestored  collisions  should  speed  up  the  Monte  Carlo 
part  of  the  program  by  a  factor  of  five  on  our  present  computer, 

b)  Higher  accuracy  Monte  Carlo  calculations  of  the  collision  integral 
by  using  more  cells  in  velocity  space  and  larger  samples  of  collisions, 

c)  Modification  of  the  program  so  that  it  will  accept  differential 
cross-sections  corresponding  to  arbitrarily  given  intarmolecular 
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potentials,  d)  Coaparison  of  our  Monte  Carlo  calculations  of  the  Boltzmann 
collision  integral  for  far-froa-equilifcriu'a  conditions  with  accurate 
and  direct  numerical  quadratures. 

Sow  let  us  consider  a  whole  class  of  kinetic  theory  problems 
for  monatomic,  single  component  gases  that  can  be  studied  directly 
with  our  present  Monte  Carlo  program  and  its  auxiliary  programs  without 
major  modification.  These  problems  are  those  in  which  there  is  just 
one  independent  variable  of  position  cr  time  in  the  Boltzmann  equation. 
Typical  problems  are  those  of  translational  relaxation  (any  such 
relaxation,  not  just  that  in  the  pseudo-shock);  shock  waves  for  a 
variety  of  Mach  numbers  and  intermolecular  potentials;  and  rarefied 
gases  in  which  there  are  large  gradients  of  temperature,  velocity,  or 
pressure  (in  other  words,  problems  of  heat  and  momentum  transfer,  of 
diffusion  and  thermal  diffusion,  and  of  gase  subjected  to  strong 
perturbations  by  high-frequency  sound  or  light). 

Useful  comparisons  will  be  possible  of  Monte  Carlo  solutions 
of  any  of  these  problems  with  algebraic  and  substitute  theories 
of  shock  structure  and  with  appropriate  approximate  theories  of  other 
physical  phenomena.  The  range  of  calculations  we  have  just  described 
could  also  be  used  as  a  basis  for  describing  phenomena  in  non¬ 
equilibrium  monatomic  gases  that  contain  a  small  fraction  of  other 
species.  Understanding  of  a  number  of  physical  phenomena  in  these 
impure  monatomic  gases  then  could  be  aided  by  the  results  derived 
from  the  Monte  Carlo  calculations  for  monatomic  gases,  namely: 
rotational  and  vibraticna!  excitation;  the  initiation  of  ionization 
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and  of  chemical  reaccions;  electrostatic  fields  produced  by  a  shock 
wave  in  a  gas  already  slightly  ionized;  and  the  phenomena  of  shock 
precursors . 

There  are  also  interesting  calculations,  which  hitherto  have 
not  been  pcssiole,  which  will  require  more  extensive  but  not  basic 
modifications  of  cur  program-  One  set  of  problems  concerns  phenomena 
in  binary  and  multi-component  mixtures  of  gases.  Treatment  of  such 
problems  requires  some  modifications  of  the  algebra  and  bookkeeping 
in  the  Monte  Carlo  program  to  account  for  collisions  between  unlike 
species-  A  second  set  of  problems  which  is  closer  to  the  interest 
of  aerodynamic is ts,  involves  two  independent  variables  of  position  or 
time.  Problems  of  this  type  are:  shock-  initiation.;  the  classical, 
problem  of  transient  heat  flow,  which  would  provide  a  useful  check 
calculation;  and  hypersonic  flow  and  heat  transfer  in  the  stagnation 
point  region  and  over  slender  bodies  for  Knudsen  numbers  less  than 
one. 
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Boltzmann  function 


Average  number  of  collisions  of  one  molecule  -  7i 

Fig.  3  Translational  relaxation  of  the  Boltzmann  function. 

The  ordinate  is  the  excess  of  the  value  of  the  Boltzmann  function, 
calculated  by  the  Monte  Carlo  solution  of  the  Boltzmann  equation,  over 
its  predicted  asymptotic  value. 
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Number  of  iterations 


Fig.  4  Convergence  of  an  iterative  solution  of  the  Boltzmann 
equation  for  a  shock  wave  (Mach  number  2.5). 
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13  abstract 

A  Monte  Carlo  method  of  solving  the  fundamental  equation  of  the 
kinetic  theory  of  dilute  gases  has  been  developed  and  successfully  applied 
to  two  problems,  one  involving  translational  relaxation  of  a  spatially 
homogeneous  gas  and  the  other  a  plane  steady  shock  of  arbitrary  strength 
(shock  strength  limited  only  by  the  fineness  of  the  velocity  space  mesh). 
This  is  the  first  and  only  method  capable  of  computing  the  molecular  veloc 
ity  distribution  under  conditions  far  from  equilibrium. 

The  essence  of  the  problem  is  evaluation  of  the  non-linear  five¬ 
dimensional  collision  integral.  Straight  forward  numerical  quadrature 
would  require  about  a  year  on  the  fastest  present  day  computers.  The 
computation  time  is  reduced  to  a  practical  value,  of  the  order  of  an  hour, 
by  a  statistical  sampling  technique  closely  resembling  the  real  statistica 
collision  phenomena  in  the  gas.  (continued  on  next  page) 
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Abstract  (continued) 


Computations  to  date  have  been  restricted  to  elastic  sphere 
scattering  of  molecules  without  internal  degrees  of  freedom.  Differ¬ 
ential  scattering  cross-sections  other  than  elastic  sphere  can  be 
accommodated  in  the  computer  program  without  complications  or  computing 
time  penalty.  Introduction  of  one  or  two  internal  molecular  degrees  of 
freedom  will  increase  the  complexity  and  computing  time,  but  not  to  an 
impractical  degree. 

Several  technical  problems  had  to  be  solved  in  order  to  make  the 
method  work  properly.  The  "fairness"  of  the  random  collision-selecting 
process  had  to  be  designed  with  care  and  verified  thoroughly.  Further¬ 
more,  a  tendency  for  the  theoretically  conserved  quantities  (number  of 
molecules,  momentum  and  energy)  to  change  slightly  because  of  interpola¬ 
tion  and  quadrature  errors,  had  to  be  corrected  to  prevent  its  building 
up  significantly  over  many  time  steps  or  iterations.  In  the  case  of 
the  more  difficult  shock  wave  computation,  a  tendency  for  the  shock  front 
to  creep  out  of  the  computational  reference  frame  in  the  course  of 
successive  iterations  had  to  be  eliminated  by  choosing  density  rather 
than  distance  normal  to  the  shock  front  as  independent  variable.  Finally, 
a  verifiably  convergent  iteration  scheme  had  to  be  devised  for  succes¬ 
sively  improving  an  initial  trial  solution.  All  of  these  technical 
problems  have  been  solved.  The  initial  trial  solution  so  far  used  has 
been  the  Mott-Smith  approximation. 

We  exhibit  for  the  translational  relaxation  problem  a  graph  of  the 
temporal  behavior  of  the  Boltzmann  function  for  Mach  numbers  ranging 
from  0.5  to  6.  For  the  shock  wave  problem  we  show,  for  Mach  number  3.0, 
contour  plots  of  the  Mott-Smith  velocity  distribution  function,  and  of 
the  collision  integral  derived  from  it,  together  with  other  plots 
characterizing  the  Monte  Carlo  solution  of  the  problem. 
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